function  signal_fft = psd2(in, npts,k)
% blackman harris windowed FFT
% removes DC component


  black_window = blackmanharris(npts);
  
  signal = (in-mean(in)*ones(1,length(in)))';
 
  signal_fft= zeros(npts,1); 
  
  for i = 1:k
      start = length(signal)-(npts-1)-(i-1)*npts;
      stop = length(signal)-(i-1)*npts;
      signal_win = signal(start:stop).*black_window;
      
      signal_fft = signal_fft + abs(fft(signal_win));
  end
    
  signal_fft = signal_fft / k;
     
  signal_fft = signal_fft / npts; % normalize